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Introduct ion 

The  design  and  optimization  of  realistic  engineering 
combustion  devices  involves  the  construction  and  execution  of 
complex  mathematical  models.  These  models  will  typically  involve 
combustion  kinetics  as  well  as  transport  processes  with  the 
physics  and  chemistry  described  by  many  parameters  which  are 
imprecisely  known.  In  addition,  existing  freedom  for  choosing 
combustion  chamber  design  will  introduce  other  potentially 
controllable  parameters  into  the  model.  Therefore,  a  central 
problem  in  all  design  problems  concerns  an  understanding  of 
system  performance  with  respect  to  its  parameter  values.  Except 
for  the  simplest  models,  such  an  understanding  will  necessitate 
extensive  computer  cal cul at  ions ,  and  repeated  execution  for  each 
new  set  of  parameters  will  lead  to  prohibitive  expense.  The  goal 
of  this  program  is  to  provide  new  insights  as  to  how  to  simplify 
detailed  submodels  which  cause  the  overall  system  calculations  to 
be  prohibitively  difficult,  and  to  exercise  the  techniques  to 
develop  simplified  chemical  kinetic  models  which  provide 
sufficient  detail  for  generating  accurate  modeling  results. 


Work  Statement 

The  work  statement  for  this  program  is  as  follows: 

1.  Model  systems  will  be  studied  to  establish  the  use  of 
elementary  sensitivity  coefficients,  Green’s  function  elements 
and  derived  sensitivity  coefficients  for  lumping  purposes. 
Appropriate  numerical  procedures  will  be  employed  including 
eigenvector-eigenvalue  analysis  and  rank  reduction  of  the  approp¬ 
riate  sensitivity  matrices. 

2.  The  sensitivity  techniques  of  item  1  will  be  developed  with 
carbon  monoxide/hydrogen/oxygen  combustion  as  a  test  case  for 
systematic  model  reduction  and  lumping.  The  ensuing  lumped 
models  will  be  compared  with  those  existing  in  the  literature  for 
their  predictive  capabilities. 

3.  Lie  algebra  techniques  for  parameter  space  mapping  and 
control  will  be  developed  starting  with  temporal  systems. 
Special  attention  will  be  given  to  using  the  techniques  for 
performing  finite  excursions  through  parameter  space.  As  the 
tools  develop,  they  will  applied  to  the  lumping  consideration 
above,  as  well  as  to  design  and  control  problems  relevant  to 
combustion  systems. 

4.  Appropriate  advanced  development  are  planned  to  extend  the 
analysis  procedures  to  more  complex  combustion  chemistry  and  to 
include  both  spatial  and  temporal  calculation  comparisons  of 
lumped  and  detailed  models. 


Status  of  Research 


A  portion  of  the  program  is  devoted  to  applying  advanced, 
local  sensitivity  analysis  and  modeling  techniques  to 
comprehensive,  elementary  reaction  mechanisms  to  aid  in  the 
development  of  such  lumped  kinetic  models.  The  chemical  systems 
of  current  interest  are  the  H2/O2/N2  and  the  CO2/H2O/O2/N2 
reactions  with  and  without  the  oxides  of  nitrogen  chemistry. 
During  the  past  year,  this  work  concentrated  on  studying  existing 
lumped  kinetic  models  in  order  to  clearly  define  and  document 
their  inadequacies  and  the  reasons  for  these  deficiencies.  In 
addition,  studies  were  initiated  to  determine  the  " lumpabi 1 ity" 
(i.e.,  the  achievable  degree  of  reduction)  of  the  chemical 
systems  of  interest.  Along  these  lines,  work  is  proceeding 
towards  the  development  of  new  lumped  parameter  models.  Details 
of  this  work  are  described  briefly  below,  and  in  more  detail,  in 
the  attached  manuscripts  which  have  been  submitted  for  review. 

Existing  lumped  parameter  models  have  been  analyzed  by 
comparing  the  predictions  of  these  reduced  models  with  those  of 
validated,  elementary  reaction  models.  Furthermore,  the  original 
semi-empirical  parameters  were  re-evaluated  from  the  data 
generated  from  the  detailed  calculations.  This  is  in  contrast  to 
the  traditional  evaluation  of  these  parameters  via  fitting  the 
empirical  models  to  measured  experimental  data.  The  results  of 
these  studies  generally  showed  existing  models  to  be  valid,  but 
only  within  the  limited  ranges  of  environmental  conditions  over 
which  they  were  originally  calibrated.  Hence,  one  may  conclude 

that  it  is  crucial  that  any  newly  developed  lumped  parameter 
model  be  validated  over  the  entire  range  of  environmental 
conditions  of  interest  to  the  specific  application.  More 

importantly,  lumped  parameter  models  valid  for  one  type  of 
physical  problem  were  generally  found  not  valid  for  another 
problem.  For  example,  a  model  constructed  from  a  purely  chemical 
system  is  not  extendable  to  one  with  transport.  Although  these 
early  conclusions  might  at  first  glance  give  a  bleak  picture  to 
the  outcome  of  the  proposed  work  (in  that  any  lumped  parameter 
model  developed  will  be  problem  and  condition  specific),  the 
results  of  gradient  sensitivity  analyses  have  revealed  a  more 
encouraging  outlook. 

Gradient  sensitivity  analysis  has  been  applied  separately  to 
both  detailed  elementary  models  and  to  lumped  parameter  models 
and  simultaneously  to  both  types  of  models  in  order  to  analyze 
for  correlations  between  parameters  from  each  type  of  model.  The 
results  from  these  studies  have  rigorously  shown,  as  a  function 
of  the  environmental  conditions,  that  changes  in  the  values  of 
lumped  parameters  could  be  directly  correlated  to  changes  in  the 
underlying  structure  of  the  elementary  reaction  mechanism. 
Furthermore,  these  results  indicated  that  to  accommodate  for  the 
changes  in  the  elementary  reaction  mechanism  structure,  both  the 
lumped  model  parameters  and  t he  governing  lumping  equations  must 
be  reformulated. 


♦  t 

Examinations  of  the  Green’s  function  response  gradients  and 
the  first  order  elementary  gradients  showed  striking  similarities 
in  their  profiles.  These  similarities  occur  for  a  number  of 
reasons  and  generally  indicate  coupled  dependencies  among  the 
variables  of  the  system.  For  example,  the  Green’s  function 
coefficients  of  the  CO/H2O/O2  reaction  showed  similar 
response  surfaces  for  ail  intermediates  in  temporal  problems. 
This  particular  observed  similarity  was  a  result  of  the  steady 
state  coupling  among  all  intermediates.  From  a  model  reduction 
standpoint,  these  types  of  couplings  can  be  used  to  evaluate  the 
degree  of  system  reduction. 

Current  research  is  investigating  systematic  procedures  to 
determine  the  maximum  degree  of  reduction  in  the  system  to 
maintain  the  accuracy  of  selected  observables.  These  procedures 
include  eigenvalue  -  eigenvector  analysis  and  rank  reduction  of 
the  gradient  matrices  described  above. 

Over  the  past  year  progress  has  also  been  made  on  three 
other  projects  all  in  general  areas  relating  to  the  mathematics 
of  kinetic  system  reduction  and  model  analysis.  These  efforts 
are  of  potentially  broad  relevance  in  combustion  and  related 
kinetic  phenomena. 

T.  The  Rigorous  Criteria  for  the  Lumping  of  Kinetic  Systems 

The  lumping  of  kinetic  systems  has  been  traditionally 
approached  by  judiciously  applying  kinetic  insights  to  reduce 
kinetic  models.  Remarkable  progress  has  been  achieved  by  this 
simple  means,  but  it  is  clear  that  a  more  rigorous  and  systematic 
approach  is  necessary.  As  a  first  step  in  this  direction  we  have 
presented  the  necessary  characteristics  to  be  satisfied  by  a 
coupled  set  of  kinetic  equations  in  order  for  them  to  admit  to 
lumping  or  system  reduction.  The  criteria  extend  previous  work 
along  these  lines  which  was  restricted  to  only  linear  kinetic 
systems.  A  set  of  theorems  has  been  established  to  test  whether 
a  system  is  rigorously  lumpable  based  on  the  form  of  the  chemical 
mechanism  or  more  explicitly,  the  Jacobian  of  the  set  of  kinetic 
equations.  The  lumping  is  defined  by  achieving  a  reduced  number 
of  chemical  species  which  are  in  general  expressed  as  a  linear 
combination  of  the  original  set.  As  a  result  of  the  theorems  and 
analysis  it  is  possible  to  establish  the  form  of  this  linear 
transformation  generating  the  lumping  process.  Although 
realistic  kinetic  systems  will  not  likely  adroit  to  exact  lumping, 
the  significance  of  this  work  resides  in  the  fact  that  we  have  a 
basis  upon  which  to  build  approximate  lumping  schemes  utilizing 
knowledge  of  the  exact  criteria  for  lumpability.  Approximate 
lumping  schemes  based  on  these  notions  are  currently  being 
developed  as  a  continuation  of  this  project. 


II.  Lie  Q§0§L9i2L§  for  the  t§L§ID2l:2r  Space  Mapping  of  Kinetic 
Systems 

The  practical  understanding  of  any  kinetic  system  requires 
detailed  knowledge  about  its  parameter  space  which  is  spanned  by 
the  system  rate  constants  and  initial  chemical  species 
concentrations.  In  some  cases,  these  parameters  may  be 
controlled  in  the  laboratory  and  in  other  situations  knowledge  of 
the  kinetic  system  behavior  in  the  space  is  sought  for  a  more 
thorough  understanding  of  the  kinetic  processes.  Research  is 
continuing  on  the  establishment  of  a  set  of  Lie  generators 
capable  of  mapping  solutions  to  the  kinetic  equations  throughout 
this  parameter  space.  In  our  previous  work  we  had  already 
established  the  differential  equations  determining  the  generator 
accomplishing  this  task  and  had  illustrated  it  for  linear  kinetic 
systems.  As  a  natural  extension  of  this  work  we  have  now  treated 
nonlinear  systems  through  second  order  kinetics.  In  particular 
we  have  shown  that  it  is  possible  under  certain  conditions  to  map 
from  first  to  second  order  kinetics  and  visa  versa.  As  a  simple 
illustration  of  these  techniques,  two  coupled  kinetic  equations 
exhibiting  competing  chemical  processes  have  been  studied  and  the 
explicit  time  independent  Lie  generators  established.  The  use  of 
symbolic  computation  has  been  quite  helpful  and  will  be 
especially  important  for  more  complex  problems.  In  the  future, 
this  work  will  be  extended  to  include  the  possibility  of 
constraints  on  the  system  parameters  and  this  situation  is 
especially  important  for  control  or  optimization  applications  of 
the  technique. 


III.  The 
Operators 


Er§£li2§I  Factorizaiion  of  Kinetic  System  Kv2lyil2D 


The  critical  step  in  both  solving  a  set  of  kinetic  equations 
as  well  as  analyzing  its  parametric  properties  using  Lie  group 
techniques  involves  the  calculation  and  analysis  of  the  evolution 
operator  of  the  physical  system.  Knowledge  of  this  operator 
essentially  provides  a  complete  description  of  the  system 
solution  and  traditional  methods  based  on  direct  numerical 
integration  of  the  equations  do  not  lead  to  the  necessary 
analytical  structure,  particularly  important  when  seeking 
physical  insight.  In  this  work  we  have  established  a  new  method 
to  factorize  the  evolution  operators  into  an  infinite  product  of 
simple  evolution  operators.  The  method  uses  Lie  operator  algebra 
and  the  evolution  operators  take  on  an  exponential  form.  The 
method  has  broad  applications  including  to  the  areas  of 
sensitivity  analysis  and  the  solution  of  kinetic  equations.  A 
sequence  of  approximants  to  the  evolution  operator  may  be 
generated  and  under  certain  conditions  the  convergence  of  these 
approximants  is  remarkably  high.  This  work  considered  the 
general  formulation  of  the  scheme  as  well  as  its  convergence 
properties.  Research  is  continuing  on  further  development  of  the 
procedure  as  well  as  applications  to  nonlinear  kinetic  systems. 
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Abstract 


In  this  work  a  new  method  to  factorize  certain  evolution 
operators  into  an  infinite  product  of  simple  evolution  operators 
is  presented.  The  method  uses  Lie  operator  algebra  and  the 
evolution  operators  are  restricted  to  exponential  form.  The 
argument  of  these  forms  is  a  first  order  linear  partial 
differential  operator.  The  method  has  broad  applications 
including  to  the  areas  of  sensitivity  analysis,  the  solution  of 
ordinary  differential  equations  and  the  solution  of  Liouville's 
equation.  A  sequence  of  f -approx imants  are  generated  to  repre¬ 
sent  the  Lie  operators.  Under  certain  conditions  the  convergence 
rate  of  the  f-approximant  sequences  is  remarkably  high.  This 
work  only  presents  the  general  formulation  of  the  Bcheme  and  some 
simple  illustrative  examples.  Investigation  of  convergence  prop¬ 
erties  is  given  in  a  companion  paper. 


I .  Introduction 


A  system  with  n-degrees  of  freedom  may  be  characterized  by 
n  variables,  xx,...,xn.  The  real  Euclidean  space  of  these 
variables  is  called  the  phase  space  of  the  system.  If  any  two 
points  in  the  phase  space  are  related  by  a  unique  transformation 
whose  functional  structure  does  not  depend  on  the  location  of  the 
points,  then  one  can  define  an  evolution  operator  for  the 
system.  In  mathematical  form  this  can  be  expressed  as 

=  Q(X,)  (1.1) 

where  xi  and  jcf  represent  the  initial  and  final  point, 
respectively.  Since  any  two  points  of  an  n-dimensional  phase 
space  may  be  connected  by  a  continuous  curve,  it  is  possible  to 
use  a  tracing  parameter  which  defines  the  position  of  the  system 
on  this  curve  during  its  evolution  from  its  initial  state  xi  to 
its  final  state  xf.  This  circumstance  often  arises  where  time  is 
the  evolutionary  parameter  and  we  will  accordingly  denote  the 
parameter  as  t.  Therefore,  the  initial  and  final  states  of  the 
system  can  be  characterized  by  the  scalar  instants  of  time  tf  and 
tf.  Hence  the  evolution  operator  Q  can  be  represented  as 
Q(t*,t,)  and 

8<  =  Q(VV  *  8,  d-2) 

where  the  dot  is  used  to  symbolically  represent  the  effect  of 
Q  on  x,  which  is  generally  nonlinear  in  character.  In  many 
applications  one  can  find  practical  expressions  for  the  operator 
Q  if  tf  and  tf  are  sufficiently  close  to  each  other.  Hence, 


the  global  evolution  operator  Q(tf,t4)  may  be  factorized 
into  a  simple  sequence  of  evolutionary  steps 

Q(tf,t1)  =  Q(tf,t.)  •  Q(t,,t.-1)...Q(t1,t1)  (1.3) 

and  by  choosing  m  sufficiently  large  this  factorization  can 
characterize  the  global  evolution  of  the  system.  If  the 
simple  short  time  interval  solutions  were  exactly  calculable, 
then  the  number  of  such  evolutions,  m,  would  lose  its  importance. 
However,  in  reality,  even  the  simple  evolutions  can  often  be  only 
approximately  determined.  In  such  a  case,  the  number  of  incre¬ 
ments  m  is  significant  since  errors  can  accumulate  to  possibly 
create  numerical  instabilities.  In  addition,  the  factorization 
requires  operators  at  times  other  than  the  initial  and  final 
specified  values.  Therefore,  a  more  global  factorization  of  the 
evolution  operator  such  as  suggested  in  this  paper  would  be  more 
attractive. 

The  present  work  considers  the  factorization  of  the 
evolution  operator  into  a  sequence  of  simple  global  evolution 
operators.  The  scheme  presented  will  maintain  its  validity  only 
on  a  special  subclass  of  evolution  operators.  First,  we  restrict 
the  system  under  consideration  to  being  autonomous  such  that  the 
evolution  operator  has  the  following  simple  structure 

Q(tf,t,)  =  Q(tf-t,)  (1.4) 

In  general,  the  complete  class  of  all  evolution  operators  is  a 
union  of  non-autonomous  and  autonomous  evolutionary  operators. 

We  also  restrict  ourselves  to  autonomous  evolution  operators 
having  an  exponential  form 
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Q(t,-tt)  =  e(t<'t<)S  (1.5) 

where  S  denotes  a  time -independent  operator.  The  class  of 
exponential  operators  can  be  divided  further  into  subclasses  by 
specifying  the  nature  of  the  operator  S.  Among  these  subclasses, 
the  quantum  mechanical  exponential  operators  can  be  defined  as 
follows 

S  =  -  i  H 

where  dimensionless  quantities  are  used  and  H  denotes  the 
Hamiltonian  of  the  system  under  consideration.  Unless  H  has 
special  structure,  the  global  factorization  of  this  type  of 
evolution  operator  using  Lie  operator  algebra  is  not  possible  due 
to  the  fact  that  H  is  generally  a  second  order  partial  differ¬ 
ential  operator  which  will  not  form  a  closed  finite  group  under 
commutation  operations.  This  difficulty  can  nevertheless  be 
circumvented  for  Hamiltonian  systems  by  consideration  of  the 
Liouville  operator.  This  latter  operator  is  amongst  the  classes 
considered  below. 

An  important  class  of  evolution  operators  is  included  in 
the  following  definition 

S  =  L  =  j  =  i  VXi . V  af^  (I,6) 

where  the  dimension  or  number  of  degrees  of  freedom  of  the  system 
may  be  finite  or  infinite.  The  finite  dimensional  case  may  be 
directly  related  to  the  corresponding  Initial  value  problem 
produced  by  the  set  of  ordinary  differential  equations1, 

4  *  f  (x  ,x  , . . . ,x  ) .  Since  almost  every  partial  differential 

j  j  1  2  N 


i.*  • 


*  *  *  »  ’  »  •  '  •  ’  »  *  »  *  »  k  »  ‘  ‘  «  *  »  *  »  •  k  ■  t,  ~  m  m 


equation  with  initial  conditions  can  be  cast  into  an  infinite  set 
of  ordinary  differential  equations  through  an  appropriately 
chosen  basis  set  expansion,  we  may  consider  the  Lie  operator  in 
Eq.  (1.6)  as  capable  of  treating  a  wide  class  of  problems.  Some 
caution  is  still  required  since  the  coefficients  in  Eq.  (1.6) 
are  scalars  while  some  formal  reductions  of  partial  differential 
equations  to  ordinary  differential  equations  can  produce  matrix 
coefficients.  In  summary,  we  restrict  ourselves  to  operators 
having  the  structure  of  Eq.  (1.6)  and  of  finite  order  N. 

Lie  operators  also  arise  in  other  areas  besides  those 
mentioned  above.  For  example,  the  investigation  of  analytic 
simplectic  maps2  and  the  description  of  the  behavior  of  tra¬ 
jectories  near  a  reference  trajectory  for  a  general  Hamiltonian 
system3  are  also  other  applications.  Recent  works4”6  have  con¬ 
sidered  the  use  of  Lie  transformations  to  perform  parameter 
space  mapping  of  the  solution  of  ordinary  differential  equations. 
Other  applications  may  also  be  found. 

The  remainder  of  this  paper  is  organized  in  the  following 
fashion.  Section  II  gives  the  general  formulation  of  the  global 
factorization  for  one-dimensional  systems  followed  by  a  generali¬ 
zation  to  multi -dimensional  systems  in  Section  III.  Some 
illustrative  examples  are  treated  in  Section  IV  and  concluding 
remarks  are  given  in  Section  V. 


Factorization  Procedure  in  the  One-Dimensional  Case 


Lie  exponential  evolution  operators  defined  by  Eqs. 

(1.5)  and  (1.6)  frequently  arise  in  many  applications.  One 
application  that  was  mentioned  above  arises  in  the  treatment  of 
ordinary  differential  equations.  In  particular,  if  we  can 
evaluate  the  effect  of  the  Lie  transformation 

Q  =  etL  ;  L  =  f(jc)  •  V  (II. 1) 


on  the  position  vector  x  around  a  point  a  in  the  phase  space  of  a 
system  defined  by 

x  =  f(x)  ,  (II. 2) 


then  the  solution  to  these  equations  may  be  written  as 


X(a,t)  =  [etL  x]x=§ 


(II- 3) 


where  x,a  and  v  are  defined  in  the  following  manner 

*T  =  Ovx2’*’xJ 

§T  =  [aifa2,...,aN] 


VT  =  r  -a  -a-  _.a— i 

~  L  ax'  ax,’ *  *  * 'ax  J 

X  2  n 

This  relation  between  the  solution  of  ordinary  differential 
equations  and  Lie  transformations  may  conversely  be  used  to 
determine  the  action  of  the  operator  Q  on  the  position  vector  by 
solving  the  following  ordinary  differential  equation 

f(x.t)  =  f(()  i  { (x,0)  =  x  (II. 7) 


(11. 4) 

(11. 5) 

(11. 6) 


This  approach  to  determining  Q  is  generally  not  preferable  since 
Eq.  (I 1.7)  is  often  only  soluble  by  elaborate  numerical  tech¬ 
niques  which  will  hide  the  important  structure  of  the  desired 


transformation.  Although  the  approach  pursued  here  is  also 
approximate,  it  will  still  leave  the  structure  of  the  evolution 
operator  rather  apparent. 


In  order  to  appreciate  the  approach  taken  below. 


some  important  properties  of  Lie  transformations 
etL[f(x)g(x)l  =  IetLf (x) 1  CetLg(x)] 

etLf(x)  =  f(etLx) 


we  recall 


(11. 8) 

(11. 9) 


The  first  of  these  equations  states  that  a  Lie  transformation  on 
a  product  of  two  functions  f(x),  g(x)  can  be  factorized  to  the 
product  of  the  Lie  transformation  on  the  individual  functions. 
This  property  is  due  to  the  exponential  structure  of  the  Lie 
transformation  along  with  application  of  the  Leibnitz  rule  of 
differentiation,  and  the  relation  is  valid  provided  that  f  and  g 
are  infinitely  differentiable  functions.  The  penetration 
property  in  Eq.  (I I. 9)  also  follows  due  to  the  particular 
structure  of  the  Lie  transformation  and  the  assumed  infinitely 
differentiable  nature  of  the  function  f.  Finally,  one  additional 
well  known  property  of  Lie  transformations  concerns  the  special 
case  of  the  translation  operator 

e'i’S  f(x)  =  f(x  +  t  a)  (II. 10) 


which  followed  from  a  simple  Taylor  expansion  of  the  right  hand 
side. 

We  now  desire  to  investigate  the  factorization  of  Lie 
transformations  for  one -dimensional  systems.  Although  the  one¬ 
dimensional  nature  of  the  problem  makes  it  formally  rather 
simple,  this  case  also  provides  the  best  means  to  develop  the 


factorization  scheme  presented  here.  In  this  case  the  Lie 
transformation  can  be  written  as 

Q  =  exp[tf ( x )  |^]  (II. 11) 

where  f(x)  may  have  a  number  of  zeros  with  one  assumed  to  exist 
at  the  origin  of  the  complex  x-plane.  This  assumption  about  the 
location  of  a  zero  of  f(x)  at  the  origin  does  not  create  any  loss 
of  generality  since  a  simple  translation  can  bring  one  of  the 
zeros  of  f(x)  to  the  origin.  The  assumption  about  the  existence 
of  at  least  one  zero  of  f(x)  is  more  restrictive.  However,  in 
problems  where  f(x)  forms  the  right  hand  side  of  an  ordinary 
differential  equation,  there  will  usually  be  at  least  one 
stationary  point  for  the  solution.  Therefore,  the  assumption 
about  the  existence  of  a  zero  of  the  function  f(x)  may  be 
regarded  as  a  minor  loss  of  generality. 

We  may  now  make  the  additional  assumption  that  the  function 
f(x)  may  be  expanded  in  a  Taylor  series 

f(x)  =  fa  fjXJ  lxl<p  (11.12) 

where  the  expansion  coefficients  f^  are  taken  as  known  from  the 
definition  of  the  system.  The  expansion  above  implies  that  the 
system  is  well  characterized,  at  least  in  a  restrictive  domain 
around  the  origin  of  the  complex  x-plane.  We  seek  the  factori¬ 
zation  of  the  evolution  operator  Q  such  that  every  factor  has  an 
Independent  contribution  in  a  fashion  analogous  to  each  term  in 


the  Taylor  series  of  Eq.  (11.12).  To  this  end  we  define  the 
flexible  factorized  structure 

Q  =  exp[tf(x)|^]  e  jii  exp[aj  (t)xJ|^]  (H.13) 

where  o^(t)  are  arbitrary  at  this  point  and  yet  to  be  determined. 

Equation  (11.13)  is  the  factorization  formula  for  the  one¬ 
dimensional  case. 

For  the  one-dimensional  case  the  factorization  in  Eq.  (11.13) 
may  seem  to  be  unnecessary  due  to  the  fact  that  the  equation  x  = 
f(x)  can  be  solved  by  the  usual  techniques  of  numerical  analysis. 
However,  in  order  to  gain  insight  into  the  more  interesting 
multi-dimensional  case,  the  present  reduced  case  presents  the 
best  way  to  understand  the  theory.  Despite  the  existence  of 
some  attempts  to  factorize  Q  by  time  ordering  techniques  with 
respect  to  t,  to  our  knowledge  there  has  been  no  factorization  of 
Q  along  the  lines  presented  in  Eq.  (11.13) 

Assuming  that  (11.13)  holds  and  the  coefficients  Oj  are 
known,  it  is  a  simple  matter  to  determine  the  effect  of  the 
operator  Q  on  x.  For  this  purpose  we  can  investigate  the  indi¬ 
vidual  effects  of  the  factors  in  Eq.  (11.13) 


Q() *x  =  exp[oj (t)xJ|^]  x 


By  using  a  simple  variable  transformation 


y  =  x"(3'l) 


(11.14) 


(11.15) 


we  may  write 


q(J)x  =  exp[-( (t)|^]  y-i/(j-i)  (11.16) 


I 


and  employ  the  translation  operator  property  of  Eq.  (I I. 10)  on 
the  y-coordinate 

Q(j)x  =  [y  -  (j-l)oj(t)]"l/(j_l)  (11.17) 


or  equivalently  in  terms  of  the  x-coordinate 


Q())x  = 


[l  -  (j-l)aj(t)xJ_1]l/()  l) 


(11.18) 


In  this  formula  the  positive  branch  of  the  root  has  been  taken. 
This  is  the  fundamental  formula  of  our  factorization  and  it 
is  valid  provided  the  argument  of  the  root  appearing  in  Eq. 
(11.18)  remains  positive.  We  are  now  able  to  evaluate  the 
effects  of  the  individual  factors  in  Eq.  (11.13).  In  applica¬ 
tions  of  Eq.  (11.13)  an  approximation  to  Q  would  consist  of 
truncating  the  infinite  product  involved. 

At  this  point  we  need  to  determine  the  coefficient 
functions  Oj .  To  this  end  we  can  use  the  following  relation 


_ 

at 


E  fi 
j  =  i  i 


i  Q(0)  =  I 


(11.19) 


which  follows  from  Eqs .  (11.12)  and  (11.13).  If  we  now  write 

Q  =  Q(1)QX  =  exp[oj(t)x  |^]  Qx  (II. 20) 

we  may  arrive  at 


i 

at 


expOoittJx!^]  [f(x)|^  -  oxx  ]exp[o1(t)x  |^]  Qx 

(11.21) 


using  the  properties  in  Eqs.  (I I. 8)  and  (I I. 9).  This  result  may 
be  re-expressed  as 


AQi 

= 


[f(exp[-olX§^]  x] 

8 

- 3 - ”  -  °  i 

lexp[-axx  x 

*  57 

(11.22) 


The  following  formula 

exp[-oi(t)  x  |j  x  =  exp(-oL(t))  (11.23) 

allows  for  a  rewriting  of  Eq.  (11.22)  utilizing  the  expansion  in 


Eq.  (11.12) 


at 


6i)  +  (f2exp(-Oi(t) )  *)  + 


Q 


l 


(11.24) 

The  operator  acting  on  Qx  on  the  right  hand  side  of  Eq.  (11.24) 
is  a  power  series  in  x.  Each  of  the  termB  of  this  series  is 
independent  and  in  the  vicinity  of  the  origin  the  first  term  will 
be  dominant.  We  desire  to  make  Qx  as  slowly  varying  as  possible 
and  therefore  demand  that  the  leading  term  in  the  series  vanish 
for  this  purpose 

6r(t)  =  fi  ■,  ox( 0)  =  0  (11.25) 

The  initial  condition  haB  been  taken  as  zero  to  make  the  simple 
evolution  operator  q(x)  unitary.  Equation  (11.24)  now  has  the 


where  ;(x)  can  be  identified  from  the  remaining  series  of 
terms  in  the  brackets  of  Eq.  (11.24)  and  f(J)(0)  is  finite. 
Exactly  this  same  logic  may  be  put  forth  to  evaluate  o2(t)  by 


successively  eliminating  higher  order  powers  of  x  in  the  differ¬ 
ential  equation.  To  construct  a  general  recursion  we  assume 
knowledge  of  the  first  n  of  the  o^'s  and  write 

Q  =  f  n  Q())}  Q  (11.26) 

lj  =  i  J  n 

which  leads  to 


T?  -  *"+1  fc  •  V°>  *  1  <“-27> 

The  time -dependence  of  f(n)(x)  is  not  explicitly  shown  and  the 
function  f(n)(x)  is  regular  at  the  origin  of  the  x-plane.  We  now 
may  write 

Qn  =  Q(n+1)  Qn+1  =  exp[on+1(t)x"+»|j]  Qn+1  (11.28) 

and  obtain 


0Q 


n  +  l 


at 


*  f (n)  [exp[-on  +  1(Oxn+1|^]x]-  on  +  l 


Qn+i 


(11.29) 

by  utilizing  again  the  properties  in  Eqs.  (II. 8)  and  (II. 9). 
Employing  the  action  of  the  factorization  operator  in  (11.18) 
gives 


0Q 


n  +  l 


at 


rf  (  n )  f  X 

-  A  1 

L  (1+no  ( t)  x") l/n 

n  +  1J 

(11.30) 
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We  now  apply  logic  analogous  to  that  leading  to  Eq.  (11.25)  and 
eliminate  the  dominant  contribution  to  the  bracketed  quantity 
multiplying  the  operator  xn+1  yielding 

°n+l  =  f(0)<°)  5  °n+l(0)  =  °  (11-31) 

where  the  Initial  value  is  again  chosen  as  zero  to  make  Q^n+1^ 
unitary.  Therefore  we  conclude 


aQ 


n  +  l 


at 


=  f (n+l) (X) 


n  +  2 


a 

ax 


n  +  l 


(11.32) 


where 


f(n+l)(x)  =  z 


(n  ) 


(l+nan+1  (t)  x") 


ZTTTr]  - 


(11.33) 


This  is  a  first  order  recursion  relation  with  the  initial 
condition 

f(l)(x)  =  [f (exp(-ai(t) )  x)  exp(Oj (t) )  -  fjx]  (11.34) 


All  of  the  o-functions  can  be  evaluated  analytically  in 
principle,  however  this  is  a  tedious  task  and  the  use  of  a 
symbolic  programming  language  such  as  MACSYMA  or  REDUCE  is 
recommended.  The  first  five  of  the  o-functions  are  given  below. 


°i(t)  =  fit 
°a(t)  =  f2  8i(t) 
o3(t)  =  f3  g2(t) 


(11.35) 

(11.36) 

(11.37) 
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o4(t)  =  [f4+  g3(t)  -  g2 


f,f. 


(t) 


(11.38) 


o5(t)  = 


fc  + 


f,f 


2  A  4 


TT  *  ITT 


8i(t)  + 


f  2f  4 


2  A  3 


TT 


TT 


g3(t)  + 


*Ui 


*rr 


82(t) 


(11.39) 


where 


gn(t)  = 


l-exp(-nf it) 


nfi 


(11.40) 


We  are  now  at  a  point  to  implement  the  factorization  scheme. 
The  essential  approximation  is  to  truncate  Eq.  (11.13)  to  a 
finite  order  thereby  producing  the  following  approximant. 


*n<*»t)  =  {iU1  Q(3)}  x 


(II. *1) 


If  the  infinite  product  representation  of  Q  given  by  (11.13) 
converges,  then  the  following  result  will  hold. 


T(x,t)  =  Qx  =  exp[tf(x)§^]  x  =  Jim  *n 


(11.42) 


Since  the  action  of  Q  on  x  defines  the  fundamental  operations 
of  concern,  we  now  focus  our  attention  on  the  ^-approx imants . 

A  recursion  relation  for  these  approximants  can  be  obtained  by 
first  noting  that 


=  (  fi  QJ1  exp[on,,(t)  xn+ii_l  x 


(11.43) 


An  application  of  Eq.  (11.18)  yields 


(11.44) 


Since  a  product  of  Lie  transformations  is  again  a  Lie 
transformation,  we  may  use  the  property  in  Eq.  (11.9)  along  with 
Eq.  (11.41)  to  conclude  that 

fn(*»t) 

£n+i(x»t)  =  r  —  u/n  (IT  45} 

[l  -  n  an+1(t)  ?R(t)J  (II. 4b) 

This  is  a  rather  simple  first  order  recursion  (difference 
equation)  whose  initial  member  is  evaluated  as  follows 


*i(x»t)  =  exp [ox (t)  x  Ij]  x  ■  x  exp(ox(t))  =  x  exp(fxt) 

(11.46) 

Although  this  is  a  simple  recursion  relation,  it  is  not  typically 
suitable  for  numerical  purposes.  Numerical  instabilities  will 
occur  if  fx  is  negative  resulting  in  excessively  small  quantities 
for  large  times  t  or  also  under  the  conditions  that  x  tends  to 
zero.  In  these  cases,  error  accumulations  may  occur  due  to  the 
truncated  arithmetic  on  the  computer.  To  prevent  this  error  we 
may  renormalize  the  ?-approximants  and  define  a  new  recursion 
relation 


*n+i 


•  n 


[l  ~  °n+ix  €R]n 


*1  =  1 


(11.47) 


where 


8n+i  =  non+i  exp(nfxt) 


(11.48) 


The  relation  between  the  new  approximants  and  the  previous  ones 
is 

?n(x,t)  =  x  exp(fxt)  (11.49) 

which  also  implies 

£(x,t)  =  exp[tf(x)|^]  x  =  {(x,t)  x  exp(fxt)  (11.50) 

Since  the  term  x  exp(fxt)  characterizes  the  linear  response  of 
the  system,  we  can  consider  £(x,t)  as  a  function  measuring  the 
deviations  of  the  system  from  its  linear  response.  We  will 
accordingly  refer  to  f  as  a  "deviation  function".  As  can  be 
easily  seen,  the  £  and  £ -approximants  have  branch  points  which 
move  on  trajectories  in  the  x-plane.  The  location  of  these 
trajectories  determines  the  convergence  regions  of  the  approxi¬ 
mants.  We  shall  leave  the  discussion  of  this  issue  and  a 
comparison  of  the  € -approximants  with  Pade  approximants  to  a 
companion  paper. 


III.  GENERALIZATION  OF  THE  FACTORIZATION  SCHEME  TO 

THE  MULTIDIMENSIONAL  CASE 

The  logic  put  forth  in  section  II  for  a  systematic  factori¬ 
zation  of  one  dimensional  evolution  operators  may  now  be 
generalized  to  multidimensional  cases.  In  this  situation,  the 
evolution  operator  acting  in  a  space  of  dimension  n  has  the  form 

Q  =  exp( t  f (x) -V)  (III.l) 

where 

xT  =  Cxlt . . . ,xN]  (III. 2) 

VT  =  ft-  t— 1 

~  [1x7' * • • '3x^j  (III. 3) 

fT  =  tfi(x) . fN(*)l  (III.  4) 

M  fW  W 

The  function  f(x)  is  assumed  to  have  a  zero  at  the  origin 

iiuo  £<i>  -  0 
*** 

and  it  is  also  assumed  to  be  expandable  in  a  multidimensional 
Taylor  series  in  the  variable  xlt...xN.  This  latter  expansion 
can  be  written  in  tensor  form  as 

fx  =  <1>f1j  Xj  +  < 2  >  f 1 j  R  XjXk  +  <3)f1Jk1  XjXkX,  +  ... 

(III. 6) 

where  the  convention  of  the  explicit  summation  over  repeated 
indices  is  used  for  convenience. 

In  the  one-dimensional  case  the  operator  exp(o1x^)  played  a 
fundamental  role  in  the  first  step  of  establishing  a  recursion 
relation  for  the  approximants .  The  same  situation  occurs  again 


here  and  we  shall  denote  this  first  degree  operator  QL  as 
taking  on  the  following  form 


Ql  =  exp(xToU)v)  (III. 7) 

<M  M  M 

where  o(0  is  a  square  matrix  or  equivalently  a  second  degree 
tensor.  The  effect  of  this  operator  on  the  position  vector  x 

mm 

is 

QLx  =  exp(oT ( 1 ) )  x  (III. 8) 

•W  mm  mm 

Since  QiX  must  be  the  system  linear  response  we  can  conclude  that 

Ojk(l)  =  t  U)fhJ  j.k  =  l . H  (III. 9) 

Henceforth,  we  shall  denote  the  linear  response  of  the  system 
evolution  by  S, 

S  =  exp(tU)f)  (III. 10) 

Using  the  definition  of  the  scalar  product  of  two  tensors  of  the 
same  order, 

A,,  4  B, .  4  =  A  o  B  (III. 11) 

Ji**,JnJi***Jn  x 

we  can  write  the  evolution  operator  in  Eq.  (III.l)  as  the 

infinite  order  factorized  product 

OO 

Q  *  Qu  n  exp[o(*)©  l(*0]  (III. 12) 

k  =  2 


*  ^  m  ^ »  1 «*  *»"_  »*  4^  S  '«•  "  *«•  ■  ■  *  m  ”  «  •  *  *  h  ^  •  ' 


■  ."V  l  .  L 


A1  V  \  V 


*V*"\.  *\*  rm^rj  A V,  .*r> 
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where  o(*)  is  a  k-th  order  tensor  to  be  be  determined  and 
l(r)  is  a  tensor  valued  operator 


OO 

'J  i) 


XJ  2  » 


J  K 


-  X<  X  4 
3  2  3  3 


(111.13) 


The  tensor  product  in  the  argument  of  each  exponential  term  in 
Eq.  (111.12)  is  itself  a  sum  of  operators  which  would  be 
difficult  to  deal  with  in  practice.  Therefore,  we  have  further 
factorized  each  individual  term  in  Eq.  (III. 12)  (except  QL) 
to  obtain 


Q=  Ql  k=2  "*  exp(°J5)a---Jk  Lli?..jk>  (III. 14) 

where  it  is  understood  that  the  coefficient  functions  o(k)  are 
now  distinct  from  the  set  in  Eq.  (III. 12).  The  starred  product  in 
this  formula  means  that  the  product  operation  is  performed  over 
the  entire  domain  of  the  j -indices.  There  iB  no  unique  ordering 
to  the  factorization  in  Eq.  (III. 14)  for  a  multi-dimensional 
case.  However,  if  we  define  the  following  operators 


n ^  =  exp[o(n)  ©  L<")] 


(III. 15) 


q(") 


1*2 


.  .  .  j 


n 


(III. 16) 


one  can  prove  that 

fq(">x}  -  {q(">x}  .  ojx2"-1] 


(III. 17) 
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Therefore,  within  this  level  of  approximation  the  expressions 
in  Eqs.  (III. 12)  and  (III. 14)  may  be  considered  equivalent.  The 
form  given  by  Eq.  (III. 14)  is  practical  since  each  of  the 
sequence  of  evolution  operators  acts  on  a  particular  coordinate 
and  degree  of  freedom. 

The  procedure  for  determining  the  o  -  tensor  is  the  same 
as  in  the  previous  section,  however  all  scalars,  (except  time) 
must  be  replaced  with  tensor  quantities  and  the  conventional 
algebra  must  be  replaced  with  tensor  algebra.  The  details  of 
these  operations  will  not  be  dealt  with  further  here,  but  the 
use  of  symbolic  programming  languages  would  be  most  helpful 
in  practice.  The  second  degree  a  -  tensor  is  given  below 
as  an  example. 

t 

(III. 18) 

The  evaluation  of  the  £  -  approx imants  can  again  be  accom¬ 
plished  by  using  the  consecutive  effects  of  the  individual 
factors  of  the  evolution  operator.  Symbolic  programming 
techniques  would  likely  be  the  best  procedure  for  determining 
the  xa,...,xN  and  t  dependence  of  the  £  -  approx imants . 


IV.  ILLUSTRATIVE  APPLICATIONS 

In  this  section,  five  problems  are  considered,  each  of 
which  exhibits  different  types  of  behavior.  For  the  sake 
of  comparison  with  the  techniques  introduced  in  the  previous 
sections,  we  have  chosen  analytically  soluble  problems  as 
described  below. 

i)  f(x)  =  1  -  e*  (IV. 1) 

From  traditional  linear  stability  analysis  arguments  this 
system  is  stable  for  x>o  and  unstable  for  x<o.  There  is  also 
only  one  steady  state  point  located  at  the  origin.  An  analytic 
expression  for  the  effect  of  the  Lie  transformation  on  x  can  be 
written  as 

£(x,t)  =  exp[tf(x)|^]  x  *  -  in  [i-( i-e~x)  e-t]  (TV. 2) 

A  careful  examination  of  the  structure  of  ?(x,t)  reveals  that  its 
branch  point  traverses  the  path  from  -»  to  +«»»  along  the  horizon¬ 
tal  axes  Tin  as  time  evolves.  Pigure  la  plots  the  exact 
deviation  function  £(x,t)  and  its  first  five  approximants  £n(x,t) 
as  defined  in  Eqs.  (11.50)  and  (11.49)  respectively  for  the 
case  x=o.i.  It  is  apparent  that  the  approximants  uniformly 
converge  to  the  true  deviation  function  as  n  increases.  The 
error  between  the  true  deviation  function  and  the  n=6 
approximant  is  shown  in  Figure  lb  where  it  is  apparent  that 
the  error  decreases  monotonically  to  an  asymptotic  value 
for  large  times.  A  similar  pair  of  plots  is  shown  in  Figure  2 


for  x=s.o.  At  this  larger  value  of  x  qualitatively  similar 
behavior  occurs  but  the  rate  of  convergence  of  the  approximants 
is  slower  and  the  peak  in  the  error  function  may  be  a  signal  of 
the  loss  of  global  convergence.  The  situation  for  negative  values 
of  x  is  different  as  shown  in  Figure  3.  Figure  3  presents  the 
case  for  x=-i.o.  The  approximants  in  this  case  seem  to  show 
oscillatory  nonmonotonic  behavior  with  regard  to  true  deviation 
function.  The  error  of  each  of  the  approximants  is  qualitatively 
similar  to  that  of  Figures  1  and  2.  At  a  sufficiently  large 
negative  value  of  x  singular  behavior  shows  up  resulting  in 
apparent  non-convergence . 

i i )  f (x)  =  l  -  e~x  (IV. 3) 

This  system  is  unstable  for  positive  x  values  due  to  the 
first  Taylor  coefficient  being  positive.  It  has  only  one  steady 
state  point  located  at  the  origin  of  the  x-plane.  The  analytic 
expression  of  the  Lie  transformation  effect  on  x  is 

J  ( x , t )  =  in{i+(e*-i)et}  (IV. 4) 

The  branch  point  trajectory  of  this  system  matches  with  the 
negative  portion  of  the  real  axis  of  the  x-plane.  The  branch 
point  moves  on  this  line  towards  the  origin  as  time  evolves 
and  reaches  there  in  the  limit  that  t-»°o  .  Figure  4  shows 
the  deviation  approximants  and  the  error  of  the  fifth  member 
for  x=o.i.  Apparent  convergence  failure  is  observed,  however 
during  a  finite  time  interval  starting  from  t*o  there  is 
temporary  convergence. 
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in)  f  (x)  =  sin  x  (IV.  5) 

This  system  is  unstable  around  x=o  for  x>ot  however  it  has 
infinitely  many  steady  state  points  and  they  alternatively  make 
the  system  either  stable  or  unstable.  Figures  5  depicts  the 
approximant  behavior  for  the  case  x=i.o.  There  is  apparent  con¬ 
vergence  behavior  in  the  figure,  however  a  peak  in  the  error 
function  may  again  be  a  signal  of  the  loss  of  global  convergence. 
It  is  difficult  to  prove  this  point  from  only  a  finite  number  of 
approx imants .  An  additional  calculation  is  shown  in  Figure  6  for 
x  *  5.0  which  is  beyond  the  second  stationary  point  of  sinx. 

This  well  behaved  nature  of  the  approximants  iB  probably  due  to 
the  fact  that  all  the  branch  points  of  this  system  are  purely 
imaginary. 


1  v) 


Stakgold  problem7 


This  problem  is  associated  with  the  consideration  of  two 
coupled  nonlinear  differential  equations  with  system  coefficients 
given  by 

fi(xi»*2)  =  x*i  "  *2  ~  *i(*?+x2) 

f2(*i**2)  =  ^*2  ~  *i  "  xa(x?+x2)  (IV. 6) 


The  analytic  expression  for  the  effect  of  the  Lie  transformation 
on  the  position  vector  is 

?i(t)  =  (x1cos  t  -  x2sin  t)  e  lxlt  7?(x,t) 


f2(t)  =  (x1sin  t  +  XjCob  t)  e  ,Xlt  r)(x,t) 


(IV. 7) 


where  X  is  assumed  to  be  negative  and  17  is  defined  as  follows 


i7(x,t) 


xi  +  x\ 

fx  -  e  -2lxltl 

1  X  1 

H  e  J 

(IV. 8) 


This  system  is  stable  as  long  as  x  remains  negative.  In  the 
case  of  positive  X  the  same  condition  again  holds  but  the 
system  does  not  have  a  steady  state  solution  and  a  limit  cycle 
appears . 

In  applying  the  method  of  Section  III  to  Eq.  (IV. 6)  we  will 
find  that  the  system  has  only  the  following  non  zero  tensor 
coefficients 


U)f 


xi 


X 


(l) 


fax 


1 


(IV.  9) 


(3) 


Xlll  = 


(3) 


1123 


_  (3) 


22X1  ~ 


_  (3) 


2222 


=  -1 


(IV. 10) 


Accordingly,  the  linear  response  term  would  be  expressed  by 
the  tensor 


S  =  exp(tU)f)  (IV. 11) 

and  elements  of  this  matrix  and  its  inverse  are  given  by  the 
following  expressions 


11  = 

cos  t 

Si  2  = 

-ext  sin  t 

21  = 

ext  sin  t 

S2  2  = 

cos  t 

(IV. 12) 

H 

1 

•w-r* 

=  e-xt  cos 

t 

s(_1) 

**12 

=  e”^^  sin 

t 

(-D 

2  1 

=  e-xt  sin 

t 

g  (  ~  3  ) 
£>22 

=  e“^t  cos 

t 

(IV. 13) 

Since  the  second  degree  Taylor  expansion  coefficients  are  zero 
it  follows  that 


o(2)  =  o 


(IV. 14) 


The  third  order  terms  are  non  zero  and  o(3)  may  be  shown  as 

t 

q(-l)  q(-l)  e(-D 

tm1m2m3m4  ®m2  l  a  »m4 1  4  at 


O) 


°1l'213,4  ■ 


S  <3>f 
s,lml  fr 


(IV. 15) 

where  the  explicit  summation  rule  over  repeated  indices  is 
employed.  After  some  tedious  algebra,  one  can  show  that  all 
elements  of  the  o( 3 ) -tensor  vanish  except  for  the  following 
four  members 


°1111  (t)  =  °1122  (t)  ~  °2211  (t)  = 

—  2  Xt 

=  °2222  (t)  =  ~~~2\ - -  “  °a(t) 


(IV. 16) 


This  result  immediately  yields  the  tensor  product 

o(3)  ©  L<3)  =  CJ3(t)  |xi  +  X3j  |xx  +  Xa  j 

(IV. 17) 

As  can  be  easily  observed  the  operators  o(3)  ©  l(3)  and 

M  M 

f ( 3 )  ©  l ( 3 )  commute  and  therefore  there  will  be  no  contribution 
from  higher  degree  terms  of  the  remainder  during  the  elimination 
of  the  operator  o(3)  0  L C 3 )  from  the  structure  of  L.  In 

M  M 

addition,  there  are  no  higher  order  terms  than  these  already 
coming  from  the  structure  of  L  itself.  Hence,  we  may  conclude 
that  the  factorization  exactly  truncates  at  its  second  order 
terms  if  we  retain  a(3)  ©  l(3)  as  a  global  second  degree 

m* 

Lie  operator.  Indeed,  if  we  write 


Qx  =  exp[xT  f T  (  3 )  •  v]  exp[o(3)  ©  l(3)]  x 

=  exp [xT  fT(i).  v]  , - - - 

lx  -  20, (t) r2 J  * 


cob  ei 
sin  e 


(IV. 18) 


then  it  follows  that 


Qx  = 


ext  1 

cos  t  -sin  t' 

xi 

1+ ( e 3 X t_  j  j  x?  +  x2 

H 

K  X 

7in  t  cos  t 

x2 

(IV. 19) 


and  the  exact  result  is  obtained.  In  this  result  we  have 
first  used  polar  coordinates 


r  =  (*i2+  *2]^  008  ®  =  XifxJ+x*] 


o(3)  ©  L ( 3 )  =  o3(t)r3 


(IV. 20) 


and  then  returned  to  the  cartesian  representation  in  Eq.  (IV. 19) 


The  result  in  Eq.  (IV. 19)  is  just  a  confirmation  of  the 
operator  algebra  introduced  earlier  in  the  paper.  We  may  still 
go  a  step  further  and  factorize  the  evolution  operator  involving 
o(3)  ©  L(3)  to  obtain 

q(2)  =  exp [o3  (t)  x 3  1^-]  exp [o3  (t)  x2^ 
exp  [a3  (t)  x3  xx  1^-]  exp  [a  3  (t)  x3 

(IV. 21) 

This  different  factorization  creates  an  error  which  is  of  the 
order  of  magnitude  of  fifth  degree  terms.  In  this  case, 
an  infinite  product  appears  which  converges  about  the  origin. 

v)  space  extension 

Consider  the  system  defined  by  the  function 

f(x)  =  ✓i-x'2“  (IV. 22) 

This  function  has  two  zeros  located  at  the  points  x=+i  and 
x=-i,  however  it  does  not  fulfill  the  requirements  for  our 
method.  In  particular,  it  cannot  be  expanded  in  a  Taylor  series 
around  these  points.  Nevertheless,  the  problem  may  still  be 
approached  by  extending  the  space  to  two  dimensions  through  the 

introduction  of  a  new  variable  in  addition  to  x  as  follows 

y  =  /x-x*  (IV. 23) 

We  can  now  define  a  new  system  with  the  descriptive  functions 


*i(x»y)  =  y 
f2(*»y)  = 


(IV. 24) 


This  new  system  satisfies  all  of  the  necessary  conditions  for 
factorization.  Therefore,  in  cases  such  as  these,  the  technique 
of  space  extension  may  make  it  possible  to  factorize  Lie  trans¬ 
formations  which  otherwise  might  not  admit  to  direct  treatment. 


V.  CONCLUDING  REMARKS 


The  basic  thrust  of  this  paper  is  the  development  of  a  new 
sequence  of  approximants  appropriate  for  time  evolution 
operators  with  Lie  generator  arguments.  Although  we  have  not 
given  convergence  theorems  for  the  {-approximants  sequences, 
the  results  in  Section  IV  are  encouraging.  Rapid  and  highly 
accurate  convergence  seems  to  be  obtained  at  least  in  a 
sufficiently  closed  vicinity  to  the  origin.  The  next  step  in 
this  work,  examined  in  a  companion  paper,  is  the  investigation 
of  the  {-approximant  singularities  and  some  convergence  theorems. 

Actual  implementation  of  the  factorization  scheme,  especially 
for  multidimensional  cases  can  involve  a  considerable  amount  of 
algebra.  The  use  of  symbolic  programing  on  the  computer  would 
likely  be  a  necessity  in  these  cases,  and  this  issue  also 
needs  further  investigation  for  its  practical  implementation.  A 
number  of  applications  of  the  factorization  may  be  envisioned  as 
suggested  in  the  introduction.  Evolution  operators  of  the  type 
studied  in  this  paper  occur  in  a  wide  variety  of  problems,  but 
perhaps  the  most  obvious  and  simple  application  would  be  to  the 
solution  of  ordinary  differential  equations.  The  possible 
attraction  here  follows  from  the  fact  that  the  approximants 
provide  a  global  solution  in  time  rather  than  the  usual 
sequential  time  stepping  procedures.  A  number  of  numerical 
issues  need  to  be  addressed  for  this  case  as  well  as  other 
applications  before  the  optimal  realm  of  utility  of  the 
scheme  may  be  established. 
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Captions 

1  Plot  of  the  exact  deviation  function  £(x,t)  and 

its  first  five  approximants  €n(x,t),  n=i,...6 
for  the  characteristic  function  in  Eq.(IV.l)  where 
x=o.i  In  figure  a,  the  last  three  approximants  are 
indistinguishable  from  the  exact  result.  Figure  b 
shows  the  deviation  function  for  the  approximant 
n=5 .  The  same  line  masks  in  figure  a  will  be  used 
in  the  remaining  {-approximant  plots. 

2.  The  same  as  figure  1,  except  x»s.o.  These  results 
are  qualitatively  similar  to  those  of  figure  1 
except  now  the  convergence  rate  is  slower  and  there 
is  a  peak  in  the  error  function. 

3.  The  same  as  figure  1,  except  now  x  =  -i.o. 

Apparent  oscillatory  nonmonotonic  behavior  is 
exhibited  with  respect  to  the  true  deviation 
function  in  figure  a. 

4.  Figure  a  exhibits  the  exact  deviation  function 
{(x,t)  and  the  first  five  approximants 

ln(x,t),  n=i,...5  corresponding  to  the  fundamental 
function  in  Eq.  (IV. 3)  at  x  =  o.x.  Figure  b  shows 
the  error  function  for  n=s  approximant.  Apparent 
divergence  behavior  is  observed  at  long  time; 
however,  during  a  finite  interval  around  the  origin 
there  is  temporary  convergence.  Here,  only  the 
fifth  approximant  goes  to  infinity.  However,  the 
first  four  approximants  also  have  branch  points 


Figure  5. 


close  to  zero  but  they  do  not  make  the  denominator 
in  the  corresponding  transformation  from  the  pre¬ 
vious  approximant  zero  when  t<io.  Some  of  the 
branch  points  are  not  even  on  the  positive  real 
axis . 

The  exact  deviation  function  £(x,t)  and  its  first 
five  approx imants  {^(x.t),  n=if...s  for  the 
characteristic  function  in  Eq.  (IV. 5)  at  x=i.o. 

The  first  and  second  approx imants  coincide  as  well 
as  the  third  and  fourth  approximants .  There  is 
apparent  convergence  behavior  in  Figure  a, 
however  a  peak  in  the  error  function  in  Figure  b 
may  signal  a  loss  of  global  convergence. 


Figure  6. 


The  same  as  Figure  5f  except  that  now  x=s.o. 
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COMPLICATIONS  OF  ONE-STEP  KINETICS  FOR 


MOIST  CO  OXIDATION 


Abstract 


The  one-step  reaction  mechanism,  CO  +  j  0^  -*  CO^  with 

d[C0]/dt  =  -kov[C0]a[H20]b[02]C 

which  is  frequently  used  in  combustion  problems  when  simplified  chemistry  is 
necessary,  is  numerically  studied  in  order  to  (i)  define  its  limitations  (and 
therefore  usage)  and  (ii)  understand  the  chemical  and  physical  reasons  for 
these  limitations.  The  analysis  is  carried  out  with  the  aid  of  a  validated 
comprehensive,  elementary  reaction  mechanism  for  moist  CO  oxidation  and  by 
specialized  sensitivity  coefficients  which  correlate  the  parameters  of  the 
global  model  to  the  parameters  of  the  elementary  model.  The  results  confirm 
many  of  the  previous,  empiricially  derived,  literature  models  and  show  the 
overall  rate  constant,  as  a  function  of  temperature,  to  exhibit  non-Arrhenius 
kinetics  and  to  be  dependent  on  pressure  and  mixture  equivalence  ratio.  More 
importantly,  models  derived  from  temporally  reacting  systems  are  shown  to  be 
improper  for  use  in  modeling  systems  reacting  with  transport  phenomena.  The 
specialized  sensitivity  coefficients  are  used  to  explain  these  complex 
behaviors  in  the  overall  model.  For  the  temporal  system,  these  coefficients 
show  that  the  global  model  must  be  able  to  account  for  dissociation  and 
equilibration  at  high  temperatures,  for  explosion  phenomena  in  the 
intermediate  temperatures,  and  for  reaction  of  carbon  monoxide  with  both  the 
hydroxyl  radical  and  hydrcperoxy  radical  at  low  temperatures.  Lastly,  methods 
for  modifying  the  existing  model  or  for  developing  a  new  model  are  suggested. 


INTRODUCTION 


The  commonly  accepted  one-step  reaction  mechanism  for  moist  carbon 

monoxide  oxidation  consists  of  the  overall  chemical  equation,  CO  +  \ 

CO^,  with  the  reaction  rate  for  CO  disappearence  defined  as  -d[C0]/dt  = 
fl  bo 

kQv[C0]  [H2O]  [0^]  .  Typically,  the  overall  specific  rate  constant  is 

expressed  in  Arrhenius  form  as  k  =  A  exp(-E  /RT) . 

ov  ov  ov 

Since  first  postulated  in  1956^,  the  global  parameters,  &ov>  ^ov>  a> 

and  c,  of  this  model  have  been  empirically  derived  by  numerous  research  groups 

2 

from  various  experiments.  Most  recently  (1985),  Lyon  et  al.  have  evaluated 

these  parameters  from  flow  reactor  data  for  dilute  mixtures  of  0^  and  H^O  with 

trace  quantities  of  CO  reacting  in  He  around  1100  K  and  1.2  atm.  Their 

results  and  the  selected  results  of  others  for  different  environmental 

conditions  are  presented  in  Table  I.  Quite  evident  from  the  table  is  a  wide 

variation  in  parameter  values.  For  example,  compare  Lyon's  results  with 

Dryer's  results  obtained  for  similar  experimental  conditions,  but  with 

different  levels  of  initial  CO  concentration.  Or  compare  all  of  the  values 

(particularly  the  overall  activation  energies)  as  a  function  of  increasing 

expe  unental  temperature.  It  is  apparent,  a  unique  parameter  set  does  not 

exist,  especially  for  the  range  of  environmental  parameters  encountered  in 

most  combustion  problems.  Regardless  of  whether  these  variations  are 

attributed  to  experimental  problems  or  to  changes  in  the  elemental  chemical 

mechanism,  this  model  continues  to  be  used  (mainly  for  lack  of  a  better  model) 

3 

in  many  applied  combustion  problems. 

The  objectives  of  the  present  paper  are  to  explore  the  complex  character 
of  this  one-step  mechanism.  In  particular,  (i)  What  are  the  limitations  of 
the  model  and  under  what  circumstances  can  it  be  used?,  (ii)  Why  does  the 
model  fail  (e.g.,  is  the  failure  a  result  of  changes  in  the  chemical  mechanism 
structure,  or  of  couplings  between  chemistry  and  transport  phenomena)?  and 
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(iii),  How  can  the  model  be  corrected  (i.e.,  what  new  methods  exist  to  modify 

and  build  future  models)?  The  paper  is  organized  in  the  order  of  the  above 

three  topics,  and  is  summarized  with  concluding  remarks. 

Finally,  from  a  more  general  perspective  this  paper  proposes  and 

demonstrates  new  mathematical  techniques  for  the  analysis  and  construction  of 

global  models.  The  term  global,  as  used  here,  is  synonymous  with  terms  such 

as  "reduced"  and  "lumped",  which  are  expressions  commonly  used  outside  the 

4 

field  of  combustion.  Hence,  the  applicability  of  the  mathematical  theory 
discussed  here  covers  considerably  more  subject  areas  than  the  combustion 
kinetics  field  for  which  it  is  demonstrated. 

LIMITATIONS  OF  THE  GLOBAL  MODEL 
Generally,  global  combustion  models  are  constructed  using  experimental 
data  from  a  wide  range  of  sources.'*  Until  recently,  this  process,  shown 

pictorially  in  Fig.  1  as  path  "a",  has  been  limited  to  the  left  hand  side  of 
the  dashed  line.  Typically,  the  methodology  of  "a"  would  consist  of  least 
squares  fitting  an  a  priori  known  function  to  measured  experimental  data. 
Advances  in  the  development  of  elementary  reaction  mechanisms,  along  with 
increasingly  more  efficient  computer  codes,  have  made  it  now  possible  to 

generate  global  parameter  models  from  the  data  bases  of  elementary  systems.^ 
In  Fig.  1,  this  process  is  represented  by  the  clockwise  perimeter  line,  paths 

"b  &  c".  At  this  point,  it  is  important  to  recognize  that  path  "c"  can  be 

interpreted  as  an  analogue  to  path  "a",  or  as  will  be  discussed  later,  can 
involve  systematic  methods  for  developing  global  models  (i.e.,  methods  for 
determination  of  both  the  governing  equations  and  parameters). 

The  development  of  elementary  models  (path  "b")  is  a  subject  outside  the 
scope  of  this  work.^  For  the  present  discussion,  it  is  sufficient  to  only 
assume  the  existence  of  such  a  model.  Currently,  fuel  molecules  for  which  the 
detailed  combustion  chemistry  can  be  described  with  confidence  have  at  most 


three  carbon  atoms;  butane  through  octane  combustion  models  are  topics  of 

9a  b 

exploratory  research.  5  Surely,  continued  progress  in  this  area  can  be 
anticipated,  but  it  should  also  be  emphasized  here  that  complete  comprehensive 
reaction  mechanisms  are  not  a  necessary  criteria  for  the  success  of  the 
present  work. 

Most  combustion  problems  may  be  described  in  generalized  form  by  the 

vector  set  of  differential  equations, 

L(Od,a)  =  0,  (la) 

where  0d(x,t,a)  is  the  vector  of  dependent  variables  (species  concentrations, 

temperature,  etc.)  desired  by  solution  of  Eq.  (la)  and  L.  (i=l,...,N  )  is  an 

1  ot 

appropriate  differential  operator  for  the  i-th  dependent  variable.  Here,  the 

superscript  "d"  refers  to  the  elementary  (or  detailed)  model,  and  the  vector  o 

includes  all  the  input  parameters  of  the  model.  The  independent  variables  of 

* 

space  and  time  are  denoted  as  x  and  t,  respectively.  Depending  upon  the 

particular  problem,  numerous  computational  codes  are  currently  available  for 

the  solution  of  Eq.  (1),  e.g. ,  for  temporally  reacting  systems,10’*1  for 

12  11 

steady  premixed  flames,  for  unsteady  premixed  flames,  for  diffusion 
flames, etc. 

In  the  present  work,  the  recent  elementary  reaction  mechanism  of  Yetter 
b 

et  al.  ’  for  describing  carbon  monoxide  /  hydrogen  /  oxygen  kinetics  is 
employed.  This  elementary  reaction  mechanism,  consisting  of  11  chemical 


A  parallel  system  of  equations  exist  for  the  global  model, 

L(0g,§)  =  0,  (lb) 

where  the  superscript  "g"  denotes  global  and  the  vector  §  includes  all  the 
input  parameters  of  this  model.  By  definition,  and  Nq  >  N^. 
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species  and  27  reversible  elementary  reactions,  has  been  validated  over  a 

combined  temperature  range  of  800  -  2800  K,  fuel  -  air  equivalence  ratios 

between  0.001  and  6,  and  pressures  between  0.3  and  3  atmospheres.  Thus,  the 

overall  rate  constant  can  be  determined  for  any  set  of  environmental 

parameters  using  numerically  generated  data  from  the  detailed  calculations. 

For  example.  Fig.  2a  shows  the  overall  rate  constant  as  a  function  of  inverse 

temperature.  Here  the  values  of  a,  b,  and  c  were  taken  from  Dryer  and  Glassman 

as  1,  0.5,  and  0.25,  respectively.  These  parameters  have  been  obtained 

experimentally  by  several  other  research  groups  as  well  (see  Table  I),  and  in 

addition,  Fristrom  and  Westenberg^  and  Dryer ^  note  that  the  same  parameters 

are  obtained  if  a  number  of  reactions  involving  the  radical  species  H,  0,  and 

OH  are  assumed  microscopically  balanced  and  the  conversion  of  CO  to  CO^  occurs 

* 

only  by  way  of  CO  +  OH  ->■  C02  +  H.  For  example,  if  the  two  reactions,  0  + 

H20  =  20H  and  02  +  M  =  20  +  M,  are  assumed  in  partial  equilibrium,  the  above 

values  for  a,  b,  and  c  are  obtained  when  the  associated  equilibrium  equations 

are  solved  for  the  hydroxyl  radical  concentration  and  substituted  in  the  rate 

equation  for  the  reaction  between  CO  and  OH.  One  should  be  very  cautious  as 

to  the  universality  of  these  parameters  since  Dryer's  experiments  were 

conducted  at  extremely  fuel  lean  conditions.  Furthermore,  many  of  the 

elementary  reactions  involving  H,  0,  and  OH  radicals  have  been  shown 

14b 

previously  to  be  balanced  only  for  fuel  lean  conditions. 

The  results  presented  in  Fig.  2a  were  obtained  from  numerous  isothermal 


The  fact  that  this  reaction  is  generally  the  slow,  rate-determining  (and 
often  dominating)  step  of  the  elementary  reaction  mechanism  has  been  used  by 
many  to  justify  the  applicability  of  the  current  one-step  reaction  mechanism. 
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calculations  using  a  dilute  mixture  of  CO,  O^,  and  H^O  in  in  which  only  the 

reaction  temperature  was  varied.  The  value  plotted  at  a  given  reaction 

temperature  represents  a  mean  value  (  /  kQv(t)  dt  /  /  dt  )  over  a  period 

during  which  approximately  90%  of  the  deficient  species  was  consumed  (CO  for  $ 

5  1  and  0^  for  <t>  >  1).  This  interval  excluded  both  induction  and  equilibrium 

chemistry.  Typically,  the  variation  in  kQv  over  this  interval  was  less  than  ± 

50%  of  the  average  value.  Obviously,  kQv  does  not  exhibit  an  Arrhenius 

dependence  over  the  complete  temperature  range  of  770  to  2500  K.  More 

importantly,  the  inferred  results  of  Eqv  over  limited  temperature  ranges  agree 

extremely  well  with  the  empirical  values  of  Eqv  obtained  from  various 
14 

experiments.  The  predicted  values  of  k  are  also  in  good  agreement  with  the 

appropriate  empirical  values.  For  example,  at  1150  K  Dryer  and  Glassman  find 

C>  =  9.7x10**  (cm3mol  *)3^4s  *  and  the  elementary  model  predicts  k^v  = 

8.6x10^  (cm3mol  1)3^s  *  and  at  1400  K  Hottel  et  al.  find  kemp  =  3.9x10** 

ov 

(cm3mol  V'4s  1  and  the  model  predicts  k^y  =  2.9x10**  (cm^mol  1)3//4s  1. 

Similar  plots  of  kQv  for  variations  in  pressure,  equivalence  ratio,  and 
water  vapor  concentration  are  shown  in  Figs.  2b,  2c,  and  2d,  respectively. 
For  variations  in  these  physical  parameters,  the  global  model  is  expected  to 
exhibit  a  constant  overall  rate  constant.  The  figures  themselves  may  be 
examined  to  study  the  range  over  which  the  global  mechanism  applies. 

The  results  of  Fig.  2  clearly  demonstrate  that  if  the  global  parameters 
of  the  one-step  mechanism  are  evaluated  at  a  specific  operating  point  (i.e., 
specific  temperature  (T) ,  pressure  (P),  equivalence  ratio  (<t>) ,  mixture 


♦The  model  predictions  reported  here  are  based  on  mixtures  studied  by  the 


composition  (c^)),  significant  errors  can  result  in  predictions  at  operating 
points  other  than  where  the  model  was  originally  calibrated.  The  most 
significant  error  would  occur  in  a  system  with  a  large  change  in  temperature. 
Moreover,  all  the  overall  rate  constants  presented  in  Fig.  2  have  been 
obtained  for  temporally  reacting  systems.  Systems  such  as  these  have  often 
been  used  to  evaluate  the  global  parameters;  the  intent  of  many  of  these 
studies  has  been  to  develop  global  models  which  could  be  later  used  to  predict 
the  same  chemistry,  but  in  more  complex  environments  (e.g.,  premixed  and 
diffusion  flames,  internal  combustion  engines,  gas  turbines,  etc.).  The 
natural  question  at  this  point  is  thus,  what  is  the  behavior  of  kQv  in  a  more 
complex  environment? 

Figure  3  shows  the  results  of  the  overall  rate  constant  obtained  from 

two  reacting  systems,  one  of  which  is  a  laminar  premixed  flame  and  the  other, 

an  adiabatic  flow  reactor.  For  both  cases,  the  initial  mixture  composition 

was  the  same.  In  this  figure,  is  not  averaged,  but  is  plotted  as  a 

function  of  increasing  reaction  temperature.  The  values  plotted  correspond  to 

the  interval  (i.e.,  positions  in  the  flame  or  reactor)  from  5  to  95% 
>v 

consumption  of  0^.  Discrepancies  exist  at  both  low  and  high  temperatures, 
which  dramatically  show  the  influence  of  diffusion  on  the  chemistry.  As  a 
consequence,  global  models  derived  in  one  environment  may  not  be  applicable  to 
another.  To  emphasize  this  point  further,  Bilger1^  has  shown  that,  in  laminar 
diffusion  flames,  kQv  exhibits  a  negative  temperature  dependence! 


Note  that  the  complex  behavior  of  kQv  vs.  inverse  temperature  for  the 
exothermic  adiabatic  flow  reactor  is  similar  to  the  dilute  isothermal  reacting 
system  of  Fig.  2a. 


The  behavior  of  k  for  numerous  other  systems/environmental  conditions 
can  be  readily  obtained  from  detailed  calculations.  The  results  presented 
thus  far  have  convincingly  shown  that  the  current  model  is  limited  and  should 
be  applied,  using  the  information  of  Figs.  2  and  3,  with  considerable  caution. 
Further  computations  support  these  findings,  but  do  not  explain  their 
occurrence.  In  the  next  section,  specialized  correlation  coefficients  are 
derived  and  analyzed  to  understand  these  anomalies. 

MODEL  ANALYSIS 

Returning  to  Fig.  1,  sensitivity  analysis  theory  has  played  a  central 
role  (paths  "d  &  e")  in  the  development  and  analysis  of  detailed  elementary 
reaction  mechanisms . ^  The  same  theories  are  applicable  to  global  model 
development  (via  path  "f"),  yet  sensitivity  analysis  has  had  only  limited 
usage  in  this  area. 

Furthermore,  the  combined  connection,  path  "d-f",  has  never  been 
explored.  The  latter  connection  enables  the  possibility  of  relating  the 
parameters  of  global  models  to  the  parameters  of  elementary  models.  Such 
information  is  useful  for  a  better  understanding  of  existing  global  schemes 
and  to  address  the  necessary  questions  for  improving  or  developing  new  models, 
e.g.,  what  lumping  parameters  should  be  employed,  and  what  should  be  the 
functional  character  of  the  lumping  equations? 

The  role  of  elementary  sensitivity  analysis  enters  as  a  means  for 
assessing  the  importance  or  contribution  of  the  various  system  parameters  J5 
with  respect  to  objectives  of  interest.  The  coefficients, 

Sij  =  3°?<x,t,£)/3e., 

provide  a  direct  measure  of  how  the  j-th  parameter  controls  the  behavior  of 
the  i-th  dependent  variable  at  point  x  and  time  t. 

The  equation  that  the  first  order  elementary  sensitivity  coefficients 


satisfy  is 
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Zn=l  (3Va°n)(80n/3V  +  C3Li/aPJ )  =  0.  (2) 

The  solutions  may  be  conveniently  expressed  in  terms  of  the  system  Green's 

1 8 

function  defined  to  satisfy  the  following  equation, 

Cl  <aVS0®>  Gnn'  =  5i„' 

whereby  the  solution  to  Eq.  (2)  may  now  be  expressed  in  terms  of  the  Green's 

,  _  18 
function, 


S  (x,t)  =  iV  /*  2dx'  /^2  dt '  G  .  i  (x ,  t  :x'  ,t ' )  3L  ,(x',t’)/3p  (4) 

ij  n  —i  X|  in  n  j 

A  similar  equation  is  derivable  for  the  second  order  elementary  sensitivity 

coefficients,  3zO?/9f$ . 3$^,  which  is  again  solvable  by  the  Green's  function 
18 

method.  Although  elementary  sensitivity  coefficients  are  not  the  subject  of 
the  present  section,  they  are  required  for  evaluation  of  the  correlation 
coefficients,  and  by  themselves,  are  a  valuable  tool  for  model  analysis.  For 
example,  an  analysis  of  Dryer's  model  shows  that  the  most  important  parameter 
on  the  CO  concentration  profile  is  the  overall  activation  energy,  Eqv,  whereas 
for  Lyon's  model,  the  most  important  parameter  is  the  CO  concentration 
exponent ,  a . 

At  this  point,  it  should  be  evident  from  the  preceding  discussion  that 
any  dependent  variable  occur ing  in  both  the  global  and  elementary  models  can 
be  represented  by  the  two  equivalent  forms, 

0f(&)  =  oj(a). 


A  functional  relationship  between  the  vectors  of  global  and  elementary 
parameters  is  consequently  implied, 


fi  =£(«)•  (5) 

The  goal  now  is  to  analyze  the  above  equation  and  to  determine  as  much 
as  possible  about  the  relationship  between  these  two  vectors.  Note  that  the 
elementary  parameters  have  to  be  chosen  at  some  operating  point  g*  in  a 
parameter  space.  (For  the  present  analysis,  it  was  assumed  that  gB  was  the 
optimum  set  of  parameters  from  Ref.  14).  The  observable  profiles  and  the 


first  order  elementary  sensitivity  coefficients  at  the  reference  point  in  o 
parameter  space,  0^(a°)  and  ( 30^/ ) a 0 »  are  readily  evaluated,  (by  the 
parallel  equations  to  Eqs .  (2),  (3)  and  (4)),  which  furthermore  implies  that 
the  jl  vector  and  its  gradient  at  the  reference  point,  &  =  (5(a#)  and 

(3f5^/3a^)ao ,  may  be  determined. 

The  method  for  calculating  these  two  quantities  consists  of  utilizing 
Eq.  (5)  by  defining  a  least  squares  (or  alternatively,  any  other  minimizing 
functional)  in  order  to  determine  the  B  vector. 

Ri  =  ’\\At  <6) 

The  integrals  in  Eq.  (6)  are  over  the  interval  (or  feature)  of  interest,  t^  < 

t  ^  t^  and  Xj  i  x  ^  x^.  The  flexibility  in  choosing  this  interval  illustrates 

that  the  present  analysis  can  be  restricted  to  only  those  aspects  of  the 

problem  of  physical  concern.  The  summation  is  over  the  number  of  experimental 

runs,  ”m"  (in  this  case,  modeling  experiments),  for  which  the  global  model  is 

o  d 

to  be  calibrated.  In  the  present  context,  the  index  associating  0*  and  (h 
with  a  specific  run,  2,  is  omitted  for  notation  convenience.  Minimization  of 
the  residual  in  Eq.  (6)  can  be  achieved  by  differentiation  with  respect  to  the 


parameter  3  ^ , 


0  =  I®  /!2dt  Z^dx  (l-(0f/0°)]  (0f/0“)  (31nOf/31nB.), 

C1  X1  1111  l  j 

j=l,2,...,Ng.  (7) 

The  solution  of  Eq.  (7)  implies  the  relation  in  (5),  although  again  it  must  be 

emphasized  that  information  on  the  elementary  model  is  at  the  reference  point 

a0.  Once  again  differentiating  Eq.  (7)  with  respect  to  the  input  parameters 
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a^,  the  following  result  is  obtained, 


I  _.  Q.  (31nB  /31na,  )  =  L.,  , 

n=l  in  n  k  ik* 


where 


Q.  =  I®  Z*2 dt  Zx2dx 
in  2=1  ti  xi 

(y)t(l-21T)(31nO®/31nBn)(31nO|/31nBj)  +  (l-y)OilnO®/31nBj31nBn)]  , 


Lik  =  l“=1  J^dt  /^dx  (l-22f)(*)  (31nO?/31n&^)  (31nCT/31nak) , 

g  d 

and  y  =  (0*/0^).  Eq.  (8)  is  a  set  of  linear  algebraic  equations  determining 
the  elements  of  the  parameter-correlation,  sensitivity  coefficient  matrix, 
93^/3^.  Note  that  the  £  vector  should  actually  carry  an  additional  index 
associating  it  with  the  i-th  observable,  and  therefore  the  present  analysis 
may  be  applied  to  any  (or  all)  observables  of  the  global  mechanism. 

When  the  fitting  implied  by  Eq.  (6)  is  sufficiently  accurate,  the  second 
term  in  the  integrand  of  Eq.  (8)  may  be  neglected  and  Eq.  (8)  can  then  be 
identified  as  having  a  least  squares  form  for  determining  the  parameter  - 
correlation  coefficients. 

Returning  to  Figs.  2  and  3,  the  most  dramatic  deviation  in  the  overall 
rate  constant,  away  from  its  "expected"  value,  is  the  non-Arrhenius 
temperature  dependence  of  Fig.  2a,  and  in  the  present  section,  the  correlation 
coefficients  described  above  are  used  to  explain  this  observed,  complex 
behavior.  Three  limited  temperature  ranges  were  studied:  (i)  1600  K  S  T  £ 
2100  K,  (ii)  900  K  <  T  <  1000  K,  and  (iii)  650  K  <  T  <  750  K.  Since  the 
largest  variations  in  the  empirical  global  parameters  occur  in  the  overall 
activation  energy  and  the  overall  preexponential  factor  (see  Table  I),  the 
present  analysis  concentrated  on  these  two  parameters;  the  parameters  a,  b, 
and  c  were  taken  from  Dryer  and  were  not  determined  from  the  least  squares  fit 
analysis.  Thus,  the  analysis  was  simplified  and  the  methodology  stressed. 
Furthermore,  the  analysis  was  carried  out  with  the  same  mixture  (in  mole 
fractions,  X(C0)  =  0.01,  X(02)  =  0.20,  X(H20)  =  0.01,  X(N2)  =  0.78)  and 
pressure  (1  atm),  varying  only  the  initial  temperature. 

The  global  model,  unlike  the  elementary  model,  does  not  have  induction 
time  kinetics.  This  complication  is  commonly  remedied  by  introducing  an 
induction  time  also  in  global  form  with  another  set  of  global  parameters.  In 
the  present  analysis,  this  period  was  defined  by 


t  =  B[C0][H201°-5[02]°-25exp(F/RT) 

where  the  parameters  B  and  F  were  determined  by  the  least  squares  analysis, 

Eqs.  (6)  and  (7).  The  concentration  exponents  were  assumed  to  be  the  same  as 

those  for  the  oxidation  process.  Generally  speaking,  this  assumption  is  not 

true  (especially  for  the  oxidation  of  hydrocarbons),  and  was  applied  here  to 

illustrate  the  sensitivity  analysis  methodology.  Again,  in  a  more 

comprehensive  analysis,  the  orders  of  the  reaction  could  also  have  been 

determined  through  the  least  squares  analysis. 

The  global  parameters  A  and  E  were  determined  from  the  CO 

ov  ov 

concentration  profile  using  a  time  interval  during  which  90%  of  the  CO  was 

consumed.  Specifically,  and  t2  of  Eq.  (6)  were  defined  as  [C0(t^)]^  = 

0 . 95 [ CO ( 0 ) ]  and  [ CO ( 1 2 ) ] ^  =  0.05[CO(0)].  Due  to  the  absence  of  the  induction 

time  in  the  global  model,  a  time  shift  was  required  to  properly  compare  the 

two  models.  This  time  adjustment  was  made  by  shifting  the  results  of  the 

global  model  to  match  the  CO  concentrations  of  both  models  at  i  [ CO ( 0 ) ] .  The 

parameters  B  and  F  were  determined  by  defining  T  as  the  time  interval  from  t=0 

to  5%  consumption  of  the  CO  concentration. 

In  the  high  temperature  regime,  the  least  squares  fit  was  conducted  with 

data  from  six  numerical  runs  with  a  combined  total  of  over  200  points  for 

comparison.  The  resulting  least  squares  parameters  for  Aqv  and  Eqv  were 
11  3  -1  3/4  -1 

1.06x10  A(cm  mol  )  s  and  21  keal/mol,  respectively.  Results  from  three 

of  the  numerical  runs  are  shown  in  Fig.  (4).  The  comparison  between  global 

and  elementary  models  is  exceptionally  good  except  at  large  reaction  times. 

Correlation  coefficients  between  these  two  global  parameters  and  the 

rate  constants  of  the  elementary  mechanism  are  listed  in  Table  II.  Recall 

that  the  elementary  model  consisted  of  54  rate  parameters;  only  those  which 

had  a  significant  influence  on  the  global  parameter  values  (i.e.,  > 

0. Ix31n& , /31na .  |  )  are  listed  in  Table  II  (and  future  tables).  Table  II 

l  j  max 
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reveals  that  both  A  and  E  are  correlated  to  many  of  the  same  and  also  to 

several  different  elementary  reactions.  This  latter  fact  would  normally  not 

be  perceived  since  Aqv  and  Eqv  are  strongly  coupled  to  each  other  in  the 

global  model  (as  determined  by  the  elementary  sensitivities).  Specifically, 

both  A  and  E  are  highly  correlated  to  the  rate  constant  for  the  reaction 
ov  ov 

between  CO  and  hydroxyl  radical.  The  overall  activation  energy  is  also  highly 
coupled  to  the  equilibrium  constants  of  the  two  chain  branching  reactions, 

0  +  H  0  •*  2 OH  (E  =  17  kcal/mol)  and  H  +  0.  -*•  OH  +  0  (E  =  18  kcal/mol). 

L.  Ql  4  8 

It  should  therefore  not  be  surprising  that  the  value  of  E  closely  resembles 
the  activation  energies  of  the  two  reactions  above,  especially  since  the 
activation  energy  for  the  rate  constant,  ^qq+qjj,  in  this  temperature  range  is 
small  (-  6  kcal/mol) . 

In  contrast,  the  overall  frequency  factor  is  sensitive  to  the  rate 

constants  of  H  +  0^  +  M  and  C02  +  H.  The  importance  of  this  latter  reaction 

on  Aqv  may  be  responsible  for  some  of  the  discrepancy  between  the  global  model 

and  elementary  model  at  large  reaction  times.  The  one-step  mechanism  was 

originally  based  on  only  the  forward  step  of  CO  +  OH  and  consequently  no 

allowance  was  made  for  equilibration.  Table  II  clearly  shows  that  the  reverse 

step  becomes  important  at  these  high  temperatures.  As  reaction  times  increase 

and/or  the  reaction  temperature  is  raised,  this  reverse  step  will  inevitably 

20 

play  a  more  dominant  role.  Recent  global  models  have  added  a  second 
reaction,  C02  •*  CO  +  $  02>  to  account  for  this  event. 

The  importance  of  the  H  +  02  +  M  -*■  HO^  +  M  reaction  at  these 
temperatures  is  surprising,  but  may  be  related  to  the  fact  that  the  H  +  0^  ■* 
OH  +  H  reaction  is  nearly  equilibrated.  In  fact,  this  reaction  has  a  small 
net  rate  of  progress  in  the  reverse  direction  allowing  the  three  body  reaction 
to  compete  favorably. 

Note  also  from  Table  II  that  for  some  reactions  of  the  elementary  model, 
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the  resulting  correlation  with  the  two  global  parameters  affects  the  overall 
kinetics  in  an  opposite  sense.  For  example,  an  increase  in  rate  constant 

correlates  to  a  decrease  in  A  and  E  .  A  decrease  in  E  increases 
OH+CO  ov  ov  ov 

the  overall  rate  of  CO  disappearence ,  whereas  a  decrease  in  Aqv>  decreases 
this  rate.  However,  the  net  effect  as  determined  by  3 Ink  /31nk„^,_„  is  an 
increase.  Moreover,  these  results  suggest  that  there  may  be  situations  where 
both  Aqv  and  Eqv  are  highly  coupled  to  a  particular  elementary  reaction,  but 
their  net  effect  on  the  overall  kinetics  may  be  zero.  Information  such  as 
this  can  be  valuable  for  guidance  towards  improving  the  present  model. 

In  Table  III,  the  induction  time  parameters,  B  and  F,  and  the  associated 
correlation  coefficients  are  reported.  The  reactions  most  correlated  to  B  are 
CO  +  OH  and  CO  +  0^  and  to  F  the  reaction  0  +  H^O.  Note  for  t,  when 
31nB/31nkj  and  31nF/31nk^  have  the  same  senses,  they  each  affect  T  in  the  same 
manner.  Consequently,  it  can  be  seen  from  Table  III  that  at  a  temperature  of 
2000  K,  the  overall  effect  of  a  perturbation  in  the  rate  constant  kH+02+M  on  T 
is  nearly  zero,  whereas  at  a  temperature  of  1600  -K,  the  net  effect  of  a 
perturbation  in  ^+q2+jij  is  to  lengthen  the  induction  time. 

The  results  of  A  ,  E  ,  31nA  /Sink.,  and  31nE  /Sink,  for  the 

ov  ov  ov  j  ov  j 

intermediate  and  low  temperature  ranges  are  tabulated  in  Tables  IV  and  V, 
respectively.  Examination  of  these  tables  immediately  reveals  a  significant 
change  in  the  elementary  structure.  These  changes  appear  in  the  global  model 
as  variations  in  the  global  parameters. 

In  the  intermediate  temperature  regime,  both  the  values  of  Aqv  and  Eqv 
increase  relative  to  their  values  in  the  high  temperature  regime.  The 
controlling  elementary  /processes  are  H  +  0^  and  H  +  0^  +  M  followed  by  CO  + 
OH.  Note  that  the  reverse  reactions  of  the  branching  reactions  are  not  as 
important  as  in  the  high  temperature  regime,  and  that  the  reaction  of  CO^  +  H 
is  not  important  at  all. 


In  the  low  temperature  regime,  Aqv  and  Eqv  decrease  relative  to  the 

intermediate  temperature  regime  values,  but  remain  larger  than  their  values  in 

the  high  temperature  regime.  In  this  temperature  range,  it  is  evident  that  H 

+  C>2  +  M  remains  competitive  with  H  +  0^  for  H-atoms,  and  CO  +  0  +  M  competes 

with  0  +  H^O  for  0-atoms.  More  importantly,  the  reaction  of  CO  with  HO^  now 

affects  A  and  E  and  dominates  over  the  importance  of  CO  +  OH. 
ov  ov 

Considering  the  entire  temperature  range,  the  behavior  of  kQv  may  now  be 

explained  by  the  results  of  the  correlation  coefficients.  Starting  in  the 

intermediate  zone,  extrapolation  of  this  rate  constant  into  the  high 

temperature  zone  overpredicts  the  actual  value.  This  apparent  decrease  results 

from  the  partial  equilibration  of  several  steps,  particularly  the  CO  +  OH  = 

CO^  +  H  reaction,  and  the  increased  importance  of  radical  -  radical  reactions. 

In  extrapolating  kQv  from  the  intermediate  regime  to  the  low  temperature 

regime,  we  find  kQv  to  be  under-predicted.  Here  the  additional  reaction  of  CO 

with  H0„  (and  with  0-atoms)  increases  k  .  The  peak  in  the  overall  activation 
2  ov 

energy  in  the  intermediate  temperature  regime  is  indicative  of  a  system 
changing  from  one  of  marginal  stability  in  the  low  temperature  regime  (high 
Eqv)  to  one  of  greater  stability  in  the  high  temperature  regime  (low  Eqv). 

GLOBAL  MODEL  DEVELOPMENT 

If  the  variations  in  the  global  parameters  are  to  be  eliminated  from  the 
current  model,  it  is  clear  that  the  functional  form  of  k  and  t  must  be 
altered.  As  mentioned  earlier,  one  advantage  of  having  available  an 
elementary  model  was  that  it  could  be  used  to  generate  the  global  model 
itself.  The  mathematical  techniques  necessary  for  such  a  process  have  been 
and  continue  to  be  topics  of  much  needed  research.  The  non-linearity  in 
combustion  kinetics  equations  make  this  a  particularly  difficult  task. 

i 

However,  several  related  techniques  currently  exist,  and  to  emphasize  their 
diversity  and  complexity,  a  few  examples  are  listed  below. 


(i)  Missing  Components  Analysis  -  In  this  analysis,  additional  parameters 
are  added  to  the  existing  functional  relationship.  The  feature  sensitivity 
coefficients  are  evaluated  with  the  newly  defined  parameters  preset  at  null 
values.  The  resulting  magnitudes  of  these  coefficients  indicate  the 
importance  of  these  new  parameters  in  the  altered  functional  relationship.  The 
method  involves  a  trial  and  error  approach,  and  although  ultimately  one  must 
actually  add  the  new  components  to  test  their  effects,  some  valuable  guidance 
can  be  obtained  from  the  missing  model  sensitivity  coefficients. 

(ii)  The  correlated  behavioral  traits  of  the  Green's  function  matrix  can 
serve  as  a  lumping  reduction  diagnostic  tool  for  choosing  lumping  parameters 
to  which  the  outputs  of  interest  respond  similarly.  The  Green's  function  of 
Eq.  (3)  is  a  generalized  memory  function  which  measures  the  stability  of  the 
i-th  observable  at  time  t  and  position  x  with  respect  to  variations  of  the 
j-th  observable  at  prior  time  t'  and  position  x’ .  For  a  temporal  system,  the 
elements  of  G  take  the  form  3(K(t)/30j  (t' ) ,  t  £  t'.  To  illustrate  the 
potential  application  of  these  gradients,  Fig.  5  is  a  schematic  diagram  of  all 
the  response  functions  obtained  from  the  elementary  model  for  a  dilute  mixture 
reacting  at  1100  K  and  1  atm.  In  the  figure,  the  complex  response  surfaces 
are  represented  by  symbols;  all  elements  represented  by  the  same  symbol  have 
similarly  behaved  response  surfaces.  Those  with  the  same  shape,  but  shaded 
in,  are  of  opposite  sign.  Blank  spaces  indicate  a  response  surface  which 
could  not  be  closely  matched  with  another.  These  results  suggest  that  one 
"effective"  radical  could  supplant  all  others  since  they  have  the  same 
response  characteristics  with  regard  to  the  major  species.  In  premixed  flame 
calculations,  the  corresponding  matrix  indicates  that  two  radicals  are 
required,  a  low  temperature  radical  and  a  high  temperature  radical. 

The  rank  of  the  Green's  function  matrix  and  its  eigenvalue  -  eigenvector 
analysis  can  be  useful  for  such  lumping.  For  example,  if  the  rank  of  £  is 
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less  than  its  dimension,  then  there  is  a  corresponding  reduced  number  of 

independent  species  implied.  In  practice,  this  issue  could  be  addressed  by  a 

rank  reduction  to  some  specified  acceptable  numerical  tolerance.  Furthermore, 

the  eigenvalues  of  £  contain  valuable  information  regarding  the  natural  time 

scales  (stability)  of  the  system.  For  lumping  concerns,  this  information  may 

be  used  to  eliminate  the  portion  of  the  system  responsible  for  short  time 

transient  behavior.  The  remaining  portion  of  the  system  would  then  be 

identified  by  the  corresponding  long  time  scale  eigenvector  of  G.  If  short 

time  behavior  was  of  concern,  the  opposite  elimination  could  be  done. 

Other  recent  mathematical  techniques  which  could  be  beneficial  to  model 

24 

lumping  include  the  Lie  algebra  approach  and  the  singular  perturbation 
procedure. ^ 

CONCLUDING  REMARKS 

This  paper  has  shown  the  commonly  used  one-step  mechanism  for  moist  CO 

i 

oxidation  to  exhibit  numerous  complexities  in  its  behavior  over  the  physical 
parameter  ranges  of  importance  to  combustion  kinetics.  Furthermore,  global 
models  found  valid  for  one  physical  system  were  found  generally  not  valid  for 

I 

another.  Correlation  coefficients,  which  relate  the  fundamental  parameters  < 

between  two  chemical  kinetic  models  (global  and  elementary),  were  derived  and 
used  to  identify  which  elementary  reactions  control  the  values  of  the  global 

I 

parameters.  This  information  was  used  to  explain  the  complex  dependence  of 
kQv  on  T.  Work  in  the  future  should  be  directed  toward  the  development  and 
application  of  systematic  mathematical  techniques,  such  as  those  introduced  in  i 

I 

the  previous  section  along  with  the  correlation  coefficients  presented  here, 
to  obtain  a  better,  universal  model  for  moist  CO  oxidation. 
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TABLE  I.  OVERALL  RATE  PARAMETERS  FOR  MOIST  CO  OXIDATION* 


Reference 

T  (K) 

A  1 
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E  2 
ov 

a 

b 
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Friedman  &  Cyphers  [1] 
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5 . 3 (9 ) 3 

20 

1 
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0 

Hottel  et  al.  [23] 
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1-2(11) 

16 

1 

0.5 

0.3 

Lavrov  [22] 
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1593 

1.8(12) 

28.3 

1 

0.5 

0.25 

Dryer  &  Glassman  [21] 

1050  - 

1200 

3.9(14) 

40 

1 

0.5 

0.25 

Lyon  et  al.  [2] 

1123  - 

1298 

1.5(10) 

23.6 

1 

0.5 

0.25 

*  -d[CO]/dt  =  Aovexp(-Eov/RT)[C0]a[H20]b[02]G 

!  ,  3  .-1.3/4  -1 

(cm  mol  )  s 
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kcal/mol 

numbers  in  parentheses  denote  powers  of  ten 
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TABLE  II.  CORRELATIONS  BETWEEN  GLOBAL  AND  ELEMENTARY  KINETIC 

PARAMETERS  AT  HIGH  TEMPERATURES  (1600  -  2100  K) 

A  =  1.06x10^  (cm^mol  s  ^ 

ov 

E  =21  kcal/mol 
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TABLE  III.  CORRELATIONS  BETWEEN  GLOBAL  AND  ELEMENTARY  KINETIC 

PARAMETERS  DURING  THE  INDUCTION  TIME  AT  HIGH  TEMPERATURES 

-35  3  1  75  7  3  1  75 

B  =  9.30x10  (cm  /molec)  '  s  =  3.83x10  (cm  /mol)  s 

F  =  7.98  kcal/mol 
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d 


System  1 
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TABLE  IV.  CORRELATIONS  BETWEEN  GLOBAL  AND  ELEMENTARY  KINETIC 
PARAMETERS  AT  INTERMEDIATE  TEMPERATURES  (900  -  1000  K) 

A 


ov 


..  in20  ,  3  -1.3/4  -1 

=  7.20x10  (cm  mol  )  s 
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TABLE  V.  CORRELATIONS  BETWEEN  GLOBAL  AND  ELEMENTARY  KINETIC 
PARAMETERS  AT  LOW  TEMPERATURES  (650  -  750  K) 
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Figure  Captions 


1.  Schematic  diagram  of  a  generalized  approach  to  global  modeling  illustrating 
the  connections  between  experiment,  detailed  kinetic  model  development, 
global  model  development,  and  sensitivity  analysis  theory. 

2.  (a)  The  overall  rate  constant,  kov,  for  a  temporal  C0/H20/02  kinetic  system 
determined  from  numerical  calculation  data.  For  all  plots,  the  units  of 
kov  are  (cm3/mol)3  '  “/s.  Here,  log  kov  =  log(d[C02] /dt)  -  log[CO]  - 
0.51og[H20]  -  0.251og[02].  To  obtain  plot,  several  isothermal  calculations 
were  run  at  different  initial  temperatures.  Note  the  non-Arrhenius 
behavior  of  kov.  At  a  temperature  of  approximately  940  K,  the  system 
crosses  the  explosion  peninsula  of  1  atm.  (b)  The  overall  rate  constant 
for  a  temporal  C0/H20/02  kinetic  system  for  various  pressures  at  1000  K  and 
1500  K.  The  decrease  in  kov  at  high  pressure  (1000  K)  indicates  the  onset 
of  slow  CO  reaction.  (c)  Variation  in  kov  for  the  isothermal  C0/H20/02 
system  for  various  equivalence  ratios  at  1500  K  and  1  atm.  For  one  curve, 
$  is  varied  by  changing  the  initial  CO  concentration.  The  decrease  in  kov 
with  decreasing  [C0(t=0)]  is  proportional  to  the  decrease  in  the  OH  radical 
concentration  overshoot.  At  very  lean  equivalence  ratios,  kov 
asymptotically  approaches  a  value  determined  by  the  equilibrium  OH 
concentration.  For  the  other  curve,  $  is  varied  by  changing  the  initial  02 
concentration.  The  decrease  in  kov  at  lean  mixtures  is  surprising  and  is  a 
good  example  where  the  sensitivity  analysis  techniques  described  could  be 
beneficial.  (d)  Variation  in  kov  for  the  temporal  C0/H20/02  system  for 
various  water  vapor  concentrations. 

3.  Comparison  of  the  overall  rate  constant  obtained  from  a  laminar  premixed 

flame  and  from  an  adiabatic  flow  reactor.  The  units  of  kov  are  (cm3/mol)3/ 
“/s.  The  initial  temperature  for  the  flame  is  298  K  and  for  the  adiabatic 
flow  reactor,  835  K.  The  latter  temperature  was  used  to  initiate  the 

reaction.  Both  systems  are  at  atmospheric  pressure.  The  initial 

composition  for  both  systems  is;  X(C0)  =  0.515,  X(02)  =  0.0987,  X(H20)  = 
0.0072,  X(H2)  =  0.0079,  and  X(N2)  =  0.371.  The  large  discrepancies  at  low 
and  high  temperatures  are  due  to  the  influence  of  diffusion  on  the 

chemistry. 

4.  Comparisons  of  the  CO  concer tration  profile  from  global  and  elementary 

models.  The  global  calculations  were  obtained  with  lumped  parameters  of 
Aov  =  1 . 064x  10 1 1  (cm3/mol)3  '  ‘‘/s  and  I'.ov  =  21.0  kcal/mol.  Initial 

conditions  of  the  three  runs  are:  P  =  1  atm,  X(C0)  =  0.01,  X(02)  =  0.20, 

X(H20)  =  0.01,  X(N2)  =  0.78,  (1)  T  =  2000  K,  (2)  T  =  1800  K,  and  (3)  T  = 

1600  K. 

5.  Schematic  diagram  of  the  moist  CO  oxidation  response  function, 
JO  (t)/80.(t’),  t  £  t'.  This  gradient  is  an  element  of  the  Green's 
function  ^matrix,  G,.(t,t')  and  is  a  generalized  memory  function  which 
measures  the  stability  of  the  i-th  observable  at  time  t  with  respect  to 
variations  of  the  j-th  observable  at  prior  time  t'.  In  the  figure, 
elements  of  similar  shape  have  similarly  behaved  response  surfaces  as  a 
function  of  t  and  t'.  Those  with  the  same  shape,  but  shaded  in,  are  of 
opposite  sign.  Blank  spaces  indicate  a  response  surface  which  could  not  be 
closely  matched  with  another.  Initial  conditions  of  the  reacting  mixture 
are:  P  =  1  atm,  T  =  1100  K,  X(C0)  =  0.002,  X(02)  =  0.028,  X(H20)  =  0.01, 
and  X(N2)  =  0.96. 
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